#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include "lauralib.h"


int main(int argc, char **argv)
{
	int i, i1, i2, n, NBINS = 20, *nb, i1_max, i2_max;
	double *t, *f, *err, *ph, P, s, r, *phb, *fb, *errb, mu, wtot, SR, SR_max, q;
	FILE *infile;

	phb = (double *)malloc(2*NBINS*sizeof(double));  	//bin phase
	fb = (double *)malloc(2*NBINS*sizeof(double));   	//average magnitude in bin
	errb = (double *)malloc(2*NBINS*sizeof(double)); 	//standard deviation of magnitudes in bin
	nb = (int *)malloc(NBINS*sizeof(int));		//number of points in bin

	infile = fopen(argv[1], "r");
	n = read3col(&t, &f, &err, infile);
	fclose(infile);

	ph = (double *)malloc(n*sizeof(double)); //phases for each data point
	for(i = 0; i < n; i++) ph[i] = 0.0;

	//subtracts off the average
	mu = getmu(f, n);
	for(i = 0; i < n; i++) f[i] -= mu;

	// folds the data
	P = 1000.0;
	for(i = 0; i < n; i++) ph[i] = getphase(t[i], 0.0, P);

	// bins the data
	for(i = 0; i < NBINS; i++) 
	{
		phb[i] = (double)i/(double)NBINS;
		fb[i] = 0.0;
		errb[i] = 0.0;
		nb[i] = 0;
	}

	// calculates average bin magnitude
	for(i = 0; i < n; i++)
	{
	//	printf("%lf\t%i\n", ph[i], (int)floor(NBINS*ph[i]));
		fb[(int)floor(NBINS*ph[i])] += f[i];
		nb[(int)floor(NBINS*ph[i])] += 1;
	}	

	for(i = 0; i < NBINS; i++)
	{
	//	printf("%lf\t", fb[i]);
		fb[i] = fb[i]/(double)nb[i];
	//	printf("%lf\n", fb[i]);
	//	printf("%d\n", nb[i]);
	}

	// calculates bin standard deviation
	for(i = 0; i < n; i++)
	{
		errb[(int)floor(NBINS*ph[i])] += pow(f[i] - fb[(int)floor(NBINS*ph[i])],2.0);  //sum of squares
	} 
	for(i = 0; i < NBINS; i++) errb[i] = sqrt(errb[i]/(nb[i] - 1.0));
//	for(i = 0; i < NBINS; i++) printf("%lf\t%lf\t%lf\n", phb[i], fb[i], errb[i]);	

// 	calculates weight "w" for each point in bin according to beginning of section 2 in Kovacs 02
	wtot = 0.0;
	for(i = 0.0; i < NBINS; i++) wtot += pow(errb[i], -2.0);
	
	for(i = 0.0; i < NBINS; i++) errb[i] = pow(errb[i], -2.0)*pow(wtot, -1.0);  //errb is now the weight of the bin (w twiddle)

	i1 = 0;
	i2 = 0;
	SR_max = -1000.0;
	q = 0.9;

	for(i1 = 0; i1 < NBINS; i1++)
	{	
		for(i2 = i1; i2 < i1 + (int)(q*NBINS); i2++)
		{
			s = 0.0;
			r = 0.0;
			
			for(i = i1; i <= i2; i++) r += errb[i];
			for(i = i1; i <= i2; i++) s += errb[i]*fb[i];

			printf("%d\t%d\t%lf\t%lf\t%lf\n", i1, i2, r, s,SR);

			SR = s*s/(r - r*r);		
			if(SR > SR_max)
			{
				SR_max = SR;
				i1_max = i1;
				i2_max = i2;
			}
		}
	}

//	for(i = 0; i < NBINS; i++) printf("errb = %lf\n", errb[i]);
	printf("%lf\t%lf\t%lf\n", phb[i1_max], phb[i2_max], SR_max);
}
